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ABSTRACT 

We present a method for beam deconvolution for cosmic microwave background (CMB) anisotropy measurements. The code takes 
as input the time-ordered data, along with the corresponding detector pointings and known beam shapes, and produces as output 
the harmonic a Ttm , a Etm , and a Bim coefficients of the observed sky. From these one can further construct temperature and Q and U 
polarisation maps. The method is applicable to absolute CMB measurements with wide sky coverage, and is independent of the 
scanning strategy. We test the code with extensive simulations, mimicking the resolution and data volume of Planck 30GHz and 
70GHz channels, but with exaggerated beam asymmetry. We apply it to multipoles up to / = 1700 and examine the results in both 
pixel space and harmonic space. We also test the method also in presence of white noise. 

Key words, methods: numerical - data analysis - cosmic microwave background 



1. Introduction 

Removal of systematic effects plays an important role in the data 
analysis of modern high-sensitivity cosmic microwave back- 
groun d (CMB) experiments, such as the WM AP and Planck mis- 
sions dJarosik et al.ll201 it iTauber et al.ll2010T) . In this work we 
concentrate on one source of systematic effects, the beam asym- 
metry. We present an efficient beam deconvolution code, called 
artDeco. It can be applied to any absolute CMB experiment 
with sufficient sky coverage. 

The required input consists of the time-ordered data (TOD), 
the corresponding detector pointings, and known beam shapes. 
The primary output of the code consists of the harmonic arim, 
aEtm, and agim coefficients of the sky. From these one can further 
construct temperature and Q and U polarisation maps. 

T his is the usual deconvoluti o n map-making problem (see, 
e.g., lArmitage & Wandeltl l2004t lArmitage-Caplan & Wandeltl 
2009 HHarrison' et al.ll201 Th . Our method differs from the earlier 
works by an efficient reformulation of the deconvolution prob- 
lem through heavy use of Wigner functions. The formulation is 
general, and does not make any assumptions on the scanning 
pattern. Similar ideas have been studied by G. Prezeau, indepen- 
dently of us (private communication). 

In the CMB literature the concept of map-making of- 
ten refers to a procedure which involves removal of corre- 
lated noise. Methods for doing this have been treated in sev- 
eral papers dKeihanen et all 12005 . 20101: iPoutanen etldlfe oOl: 



Ashd own et al] 12007 a b: Kur ki-Suonio et al . 2009; S utton et alJ 
2009). The methods discussed in these works construct pixelised 
temperature and Q, U polarisation maps, but do not attempt to 
correct for an asymmetric beam shape. 

Our method steps in at the point where the correlated noise 
has already been cleaned from the data. The input TOD is as- 
sumed to include signal and uncorrected noise only. 



Since the aim of this paper is to present a new deconvo- 
lution method, rather than assess the performance of any par- 
ticular instrument, we test the code with a somewhat idealised 
simulation. In particular, we ignore other systematic effects such 
as frequency response. On the other hand, we choose strongly 
asymmetric beam shapes, in order to show the beam effects more 
clearly. 

2. Beam deconvolution 

2. 1 . Constructing the linear system 

In this section we present the deconvolution problem using a 
formalism based on Wigner functions. We write the solution in 
a form which allows to solve the problem numerically in an effi- 
cient way. Some details are left for Appendix lAl 

We start by writing the signal tj received by a detector at a 
given time Tj as 



(1) 



slink 



The a s \ m (s = 0, +2) are harmonic coefficients which represent 
temperature and polarisation of the sky signal, b s [k are corre- 
sponding beam coefficients, rij represents noise, and D' mk (u>j) is 
a Wigner function which defines a rotation of the beam from a 
fiducial orientation at the north pole to its actual position and 
orientation. For exact definitions, see Appendix lAl 

We have adopted a compact notation where a> denotes a com- 
bination of three angles {<f>, 6, i//}, which define the detector's ori- 
entation. Angles 6 and <p define a point on the celestial sphere, 
and ip defines the beam orientation. The angles can be interpreted 
as Euler rotation angles. 

We perform a simple linear least-squares fit of a s i m to the 
data, up to some multipole / max . 
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In general, a linear least-squares fit leads to an equation of 
the form 

x = (A l -C- n 1 A + C- l y 1 A l C- n i y, (2) 

where x is the vector of unknowns, y is the data vector, C„ is 
noise covariance and C s repre sents a priori information on the 
covanance of x dWiener| [l949). Matrix A relates the unknowns 
to the data through 

y-Ax + n. (3) 
In the present case matrix A is 

A L = Tj b '* D '^ oj j ) - (4) 

k 

and x represents a s i m . 

In this paper we make two simplifying assumptions, which 
however are not of fundamental nature. Firstly, we assume the 
noise is white with constant variance, and assign equal weights 
to all samples /'. Secondly, we use no a priori information on 
the sky signal, i.e. the harmonic coefficients are assumed to be 
uncorrelated and can host a signal of infinite variance. Under 
these assumptions the covariance matrices drop out of the equa- 
tion (C„=const, Cj l = 0). The least-squares solution for a s i m is 
obtained by solving the linear system of equations 

Z ^t A ili« A il'm' as 'l' m ' ~ /Ailrrfi- ^) 

s'l'm' j j 

In the following we process Eq. (O further to bring it into 
a form where it can be solved in an efficient way. We treat the 
right- and left-hand sides separately. 

2.1.1. Right-hand side 
The right-hand side yields 

Z A ^ = Z b M D U"j)tj = Z b A- ^ 

j kj k 

Here we define 

T l mk = T J D l mk (<o j )tj. (7) 

To numerically evaluate the transform, we replace the sum over 
TOD samples by a sum over a 3-dimensional grid as follows. 
We divide the space spanned by all possible triplets {</>, 6, ijj] into 
N bins of equal volume and denote by f(o>,-), i = 1 . . .N, the 
sum of valu es f, falling into the bin. In 6, <p we use HEALPijfl 
pixelisation (Gorski et al. 2005), and in we divide the space 
uniformly. 

We call quantity t(a>) the 3-dimensional signal map, and 
quantity 

^ = Z D '^ (w)f(w) (8) 

the Wigner transform of the signal map. 

Here 2 W denotes a sum over all bins Expansion ([8} can 
be evaluated numerically quite efficiently using FFT technique 
and fast algorithms for the generation of Wigner functions. 

Quatities t(u>) and Tf, do not carry memory of in which or- 
der the sky was scanned. They just record which pixels of the 
sky were observed, in which orientation of the beam. 

1 http : //sourcef orge . net/pro j ects/healpix 



2.1.2. Left-hand side 

Consider then the matrix on the left-hand side of Eq. (O, 

Z KiAi-m- = Z ^ D >AW^r (9) 

j jkk> 

The multiplier matrix is too large to be inverted. Instead, 
we use the conjugate gradient (CG) iteration method to solve 
Eq. ©. 

We pixelise the 3-dimensional pointing space in the same 
way we did with the signal, and denote by n(cS) the number of 
hits into bin co. In line with the definition of the 3D signal map, 
we call n(a)) the 3-dimensional hit map. 

We replace the sum over samples j by a sum over bins and 
write 

Z D '^j)Dlk^j) = Z D'^Dl^Maj). (10) 

j u 

As the next step, we write the Wigner functions as 

D' mk (e, 0, «A) - e- im W mk {9)e- ik ^ (11) 

where d l , are the reduced Wigner functions. Inserting this into 
( fTOb we get 

Z D »*>^'^ 

j 

= 2 d' mk (0)d l m , e (6) e^-'^e'^nie, <p, f) 

8 <t>,tp 

= Z4 t (0)4<*<(«-w(0). d2) 

e 

where we have defined 

N mk (0) = ^ e- in "fn(0, 0, ifte-**. (13) 

Substituting this back into (O and rearranging, we bring the left- 
hand side of the deconvolution equation into the form 

Z Z A l A ^m<^« < 14 ) 
s'l'm' j 

= Z b slk Z d U8) Z N m-W(ff) Z 
k 8 m'k' 1' s' 

This can be evaluated much more quickly than the original for- 
mula (0, since the summation over the long TOD vector has 
been replaced by an operation on the more compact object. 
Also, the convolution operations in indices k, m can be carried 
out efficiently using FFT technique. 

The terms are written in the order in which they are applied 
in the code. 

2.2. Formulation through 3j coefficients and preconditioner 

Eq. (0 can be written in a yet more compact form. 

The product of two Wigner functions can be expressed in 
terms of 3 j symbols (Appendix lA.3b as 

D l ^)D' mtk ,(co) = £(2/ 2 + l)D^_ mk ,_ k {co) (15) 

h 

x(-l)'"' +A "( 1 L , h ){[ V v v h X 
I m —m m — m \k -k k - k 
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Substituting this into (TlOb and that further into (0 we arrive 
at the formula 



slm s'l'm' 

j 

=Z^«'Z (2Z2+1) ^' 

kk' 

X(-l) 



-m,k'—k 

kk' 

m'+k' III' h \ I I I' h 



(16) 



m —m' m' —m)\k —k! k' — k }' 



Indices / and I' run from to l max . The sum over h covers the 
values for which [/ — 1'\ < 1% < I + /'. We have defined the Wigner 
transform of the hit map as 



(17) 



We use the formula ( fTBI l to explicitly construct the diagonal of 
the matrix, which serves as a preconditioner to speed up the CG 
iteration. 

The formulation presented here is general, and makes no as- 
sumptions on the scanning pattern, or on the sky coverage. 

3. Implementation 

3.1. User interface 

The artDeco code takes two lists of input data objects contain- 
ing the beam coefficients b s n and the so-called "3D maps" for 
the respective detector, for all detectors that are to be processed 
simultaneously. The beam coefficients can be read from FITS 
files in the format that is also being used by the HEALPix pack- 
age for storing spherical harmonic coefficients. 

The 3D maps contain the quantities n{u>) and t(u>) for a given 
detector; they are created from signal and detector pointing time 
streams by another standalone module called todtobin. This 
module achieves the angular discretisation of the detector point- 
ings by sorting them into pixels of a HEALPix map with a user- 
defined A^side parameter, which are in turn further subdivided into 
Hi/, bins for the discretisation in tfr. The choice of A s jd e and has 
direct consequences on the run-time of the deconvolver and the 
accuracy of its results: if high values are chosen, the 3D map 
will become very large, and the deconvolver will require more 
memory and CPU time. On the other hand, too low values will 
result in noticeable discretisation errors and inaccurate output 
axim- As a guideline, A s ;d e should be chosen in such a way that 
the smallest features of all beams are well resolved at the angular 
resolution 

So « 3500'/7V s ide (18) 

associated with this A s ;d e - In a similar fashion, should be large 
enough to resolve azimuthal features of the beams (assuming 
the beams are centered on a pole). For completely axisymmetric 
beams, — 1 is sufficient. Generally, n$ = 10£ max (where k max 
is the highest azimuthal multipole used to describe the beams) is 
a good choice. 

artDeco takes several user-supplied parameters influencing 
its behaviour. Most important among these are / max , which spec- 
ifies the maximum multipole up to which deconvolution is car- 
ried out, and k m . dx < Z max , which limits the maximum azimuthal 
multipole used for the beams. The explicit introduction of k lmx 
is advantageous since most beams can be expressed accurately 
even when using a £ max <K Z max (even £ max = if the beam has ro- 
tational symmetry and is unpolarised), and the reduction in k max 
reduces CPU and memory requirements dramatically. 



artDeco's output consists of a set of a; m in HEALPix- 
compatible format. Depending on whether the user requested 
polarisation, the file contains only ajim or additional a^hn an d 
asim coefficients. 



3.2. Algorithm overview 

The deconvolution algorithm has been implemented in C++, 
making use of the MPI library for parallelisation. The quantities 
which depend on 9, such as d and N, are distributed to MPI tasks 
according to the 9 angle. This allows a fairly good load balance, 
and summation over 9, when needed, is done efficiently through 
the collective MPI_Reduce() routine. The a s i m coefficients are 
distributed according to index m. The beam coefficients b s n are 
fully present on all tasks. 

The procedure of evaluating Eq. (TT4T > can be summarised 
as follows. First we multiply the input a by the beam coeffi- 
cients b. Then we perform a loop over index m', inside which 
we broadcast the a. 5 </< m < coefficients for the current m' to all pro- 
cesses, multiply by b, construct the reduced Wigner matrices for 
the local values of 9, and accumulate over The multiplication 
by matrix A is a convolution operation in indices m, k and can 
be performed efficiently using 2-dimensional Fast Fourier trans- 
forms. Each MPI task performs the multiplication by N inde- 
pendently, for the local 9. After that, we make another loop over 
index to, construct again the Wigner matrices, sum over 9, col- 
lect the result according to to to the owner process, and multiply 
by the beam coefficients. 

The solution is assumed to have converged sufficiently as 
soon as the residual of the current a s i m estimate is smaller than 
10~ 12 times the residual of a vector containing all zeroes. 

To generate the 3 /' symbols needed to construct the precon- 
ditioner, we use a C implementation of the recursive algorithm 
described bv lSchulten & Gordonl (119751) . 

For a more detailed technical description of the algorithm, 
see AppendixlBl CPU and memory requirements of various runs 
are presented and discussed in Appendix ICl 



4. Simulation setup 

The deconvolution algorithm was tested using mock data 
produced by the Planck simulation package (Level-S, 
iReinecke et al1l2006l) . 



4.1. Mission parameters 

We assumed a scanning strategy which caused the telescope's 
spin axis to perform two cycloidal excursions from the ecliptic 
plane per year with an amplitude of 7.5°; in azimuthal direction, 
the spin axis always pointed towards the Sun. 

Time-ordered data were generated for a time span of 10000 
hours, corresponding to slightly more than two full sky surveys; 
their length is on the order of 10 9 samples. 

4.2. Detector geometry 

For the first set of experiments, data was simulated for two horns, 
i.e. two pairs of detectors, which were modeled roughly follow- 
ing Planck's 30GHz channel. In particular, the detectors had an 
average FWHM of 32'.2, the polarisation directions in each de- 
tector pair were rotated roughly (but not exactly) 90° against 
each other, and the polarisation directions of the two horns were 
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Table 1. Overview of the detector parameters 



Name 


Frequency 


ellipticity 


FWHM 






30 ,o 


30GHz 


1.7 


32.2352' 


-22.3° 


514 


30 0J 


30GHz 


1.8 


32.1377' 


67.6° 


514 


30i, 


30GHz 


1.7 


32.2352' 


22.3° 


514 


30 u 


30GHz 


1.8 


32.1377' 


-67.6° 


514 


70 ,o 


70GHz 


1.7 


13.0328' 


22.2° 




70m 


70GHz 


1.8 


12.9916' 


-67.5° 




70i, 


70GHz 


1.7 


13.0328' 


-22.2° 




70u 


70GHz 


1.8 


12.9916' 


67.5° 





Notes. The first subscript denotes the horn, the second one the num- 
ber of the detector inside the horn. cr w is the standard deviation of the 
Gaussian random numbers that are added to each TOD sample in simu- 
lations with white noise. 



in turn shifted by 45°. Both horns followed the same scan path 
on the sky, although with a little time delay. 

The beam shapes were assumed to be elliptical, with elliptic- 
ities of 1 .7 and 1.8 for each detector in a pair, respectively. These 
values are considerably higher than in a typical experiment, and 
were chosen to demonstrate the performance of the method even 
under unfavourable conditions. 

In order to show the algorithm's performance at higher reso- 
lutions, another data set was produced using two horns sensitive 
at 70GHz and with an average FWHM of 13'. 

TableQ]gives an overview over the detector properties. 

4.3. Input data 

To generate a theoretical power spectrum for the CMB emission, 
we used the CAMB0 code with the cosmological parameters 
Qi = 0.045, Q cdm = 0.255 and = 0.7. Using the HEALPix 
tool syn_alm_cxx, we created a random, unconstrained realisa- 
tion of polarised a/„, from this power spectrum. 

Emission maps of foreground components were gen- 
erated using the Plan ck pre-launch Sky modefl (PSM, 
iDelabrouille et alJ 120121) . version 1.6.6. We chose to include 
the contributions from galactic infrared emission (galactic syn- 
chrotron radiation, free-free emission, as well as thermal and 
spinning dust). 

The signal from strong point sources was taken into account 
as well; in order to achieve this, a PSM-generated catalog of 
point sources was directly expanded into spherical harmonic co- 
efficients, using a cutoff value / max high enough that the beam 
response was negligible at this scale. 

Polarised emission was disabled for all contributions except 
for the CMB, which makes it very easy to detect leakage of the 
strong temperature signal in the galactic plane into the polarisa- 
tion signal. 

We neglected bandpass effects, and assumed an identical 
delta frequency response for all detectors. 

The power spectra of the signal components at 30 GHz are 
plotted in Fig.Q] 

4.4. Noise 

In some of our simulations, detector noise was added to the ide- 
alised sky signal. Since the deconvolution method requires cor- 




1000 



Fig. 1. Spectra of the 30GHz input sky components: CMB 
(blue), diffuse foregrounds (purple), point sources (red), and to- 
tal (black). 

related noise to be removed from the time-ordered data before- 
hand, we only added Gaussian white noise with an appropriately 
chosen <x to the simulated samples (see Table[TJ. 



4.5. 3D map generation 

The mock TOD were processed with the todtobin tool to gen- 
erate the 3D input maps required by the deconvolver. We set 
7Y si de = 1024 for the 30GHz case, and 7Y side = 2048 for 70GHz; 
was always set to 256. 



5. Results 

5. 1 . Constructing a sky map 

The primary output of the deconvolution code consists of the 
o-xim (X=T,E,B) coefficients of the sky, up to / max . From these one 
can construct the usual temperature and Q,U polarisation maps 
through harmonic expansion of the coefficients. We do this with 
the alm2map_cxx tool of the HEALPix package. 

For comparison, we have also produced binned maps directly 
from the TOD. A binned map is constructed by assigning each 
TOD sample completely to the pixel containing the center of the 
beam. The binning operation can be formally written as 



(p T c- n l py l p T c- n 



y- 



(19) 



http : //camb . info 



http : //www . ape . univ-paris7 . fr/- delabrou/PSM/psm. html 



Here P is a pointing matrix, C„ is white noise covariance, and m 
is a (I,Q,U) map triplet. 

Binned maps are distorted by beam effects, if the beam 
is asymmetric. In particular, the polarisation component of a 
binned map is contaminated by temperature leakage through 
beam shape mismatch. We aim at demonstrating that deconvolu- 
tion can reduce these undesired effects. 

Since our simulation includes point sources, the input axim 
coefficients are not band-limited. However, we can only solve 
the coefficients up to some limited / max . This leads to ringing 
artifacts, if a map is constructed directly from the raw axim coef- 
ficients. Ringing is particularly prominent around point sources 
and other sharp features. 

Ringing is caused by the sharp cut-off of the spectrum at Z max , 
and it occurs even if the deconvolved coefficients are fully cor- 
rect up to / max . For this reason it is necessary to smooth the co- 
efficients with a symmetric Gaussian beam prior to constructing 




Fig. 2. Ringing due to insufficient smoothing. A map constructed 
from the input aj\ m coefficients with Z max =800 with no smooth- 
ing (left) or smoothed to FWHM=32' (right). The cut-off of the 
spectrum at Z max leads to strong ringing effects around a point 
source. 




a map. The required amount of smoothing depends on the spec- 
trum of the underlying sky, and on Z max . The highest Z max that one 
can solve in turn depends on the beam width. In other words, it 
is not possible to recover structures significantly smaller than the 
beam itself. This gives the following rule of thumb: the required 
smoothing width is given by the smallest features of the original 
beam. 

In Fig. [2] we demonstrate the ringing around a point source. 
We show an image of a point source, constructed from the in- 
put axim, but including only coefficients up to Z max =800. Ringing 
artifacts are evident. On the right we show the same region of 
the sky, after the map is smoothed with a FWHM=32' Gaussian 
symmetric beam. Smoothing removes the ringing, but naturally 
also smoothes out the smallest structures in the map. 



-0.0010 ^^^^^^^^^k t^^^^^m 0.0010 

Fig. 3. 30GHz temperature map, without noise: binned (top), and 
deconvolved and FWHM=32' smoothed (bottom). 

T binned T deconvolved. fwhm=33 




5.2. Signal-only simulations with elliptic beam 

We analyze in detail first the noise-free 30GHz simulation. A 
detailed description of the simulation was given in Section |4] 
We run the code with various combinations of parameters Z max , 
k max , and analyze the results both in pixel space and in harmonic 
space. We include data from all four detectors and consider both 
temperature and polarisation. 



-0.0010 0.0010 -0.0010 H 0.0010 

Fig. 4. A zoom into the 30GHz temperature map without noise: 
binned (left) and deconvolved (right). Shown is a 1000'xlOOO' 
patch of the sky. The deconvolved map can be compared to the 
map constructed from the input aj\ m (right panel of Fig. [2J. 



5.2.1. Temperature map 

In Fig. [3] we show a full-sky Mollweide projection of the binned 
temperature map, together with that of the deconvolved map. 
In the case shown we used deconvolution parameters / max =800, 

We have chosen the value FWHM=32' as the baseline 
smoothing width for the 30GHz maps. This corresponds to the 
mean width of the input beam. In some cases it was necessary to 
apply even stronger smoothing. 

Difference between the binned and the deconvolved map 
does not appear dramatic in the full-sky map image, but becomes 
evident when we zoom into a point source. Figure [4] shows a 
randomly selected strong point source below the galactic plane. 
While the source in the binned map is clearly elongated due to 
beam shape, it appears circular in the deconvolved map. 

The deconvolved image can be compared to the right panel 
of Fig. |2l where we have shown the same patch of the sky, as 



constructed from the input ax\ m coefficients. The maps are nearly 
identical. 

5.2.2. Polarisation maps 

The beam effects show up more dramatically in the polarisation 
maps (Fig. |5). The temperature signals recorded by two detec- 
tors differ due to the different beam shapes. The binning proce- 
dure interprets this erroneously as polarisation signal. This leads 
to a leakage of temperature signal into the polarisation maps. 
Since our simulation did not include polarised foregrounds, all 
the galactic structure seen in the binned map is temperature leak- 
age through beam mismatch. 

We show deconvolved maps with two levels of smoothing 
below the binned map. Smoothing with a FWHM=32' beam 
is not sufficient to fully remove the galactic residual, but if we 
smooth the map further to FWHM=45', the residual disappears 
almost completely. 
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-1.0e-05 ^^^^^^^h ^^^^h 1.0e-05 



Fig. 5. 30GHz Q polarisation maps, in absence of noise: 
binned (top), and deconvolved with FWHM=32' (middle), or 
FWHM=45' (bottom) smoothing. All the galactic structure is 
fake polarisation signal arising from beam mismatch. 

In Fig. [6] we show a zoom into the Q polarisation map, just 
below the galactic plane. In the binned map we see a strong 
galactic residual, and typical "four-leaf clover" structures that 
arise from temperature leakage at the location of a point source. 
We show the deconvolved map again for two levels of smooth- 
ing: FWHM=32' and FWHM=45'. The additional smoothing re- 
moves the clover structure nearly completely. In order to show 
that this is not only an effect of more aggressive smoothing, we 
smoothed the binned map further with a Gaussian beam of width 
FWHM=32', which, combined with the 32' of the original beam, 
gives a total smoothing of 45'. This map is shown below the un- 
smoothed binned map. The additional smoothing has little effect 
on the galactic residual. 

5.2.3. Harmonic space 

We proceed to analyze the deconvolution results in harmonic 
space. We compute the TT, EE, and TE spectra of the decon- 
volved axim coefficients, and compare them to the corresponding 
input spectra. The spectrum is constructed as 

' = 2TTT 2-i axlmaYlm 

m=—l 

where X, Y stand for T and E. 



Q binned Q deconvolved 




Fig. 6. A zoom into the 30GHz Q polarisation map without 
noise: binned (left) and deconvolved (right). The maps were 
smoothed to resolution FWHM=32' (top), or to 45' (bottom). We 
show a 2500'x2500' patch of the sky. The horizontal structure is 
spurious signal due to temperature leakage from the galaxy. 



We run the code for all combinations of k max =4,6 and 
Z max =400,500,600,700,800. Results for different combinations 
are shown in Fig. [7] We apply no smoothing to the coefficients. 

We observe that the optimal value for Ic m . dx is coupled to 
/max- Results for k max =4 and k m . dx -6 begin to diverge only above 
/ = 400. The results for different l m . dx with k max =6 are practi- 
cally on top of each other, except for the highest multipoles near 
^max, which are clearly distorted. When constructing a map, we 
have cut out the last 20 multipoles before making the harmonic 
expansion. 

For comparison we show also the spectrum obtained from 
harmonic expansion of the binned map, which we compute using 
the anafast tool of the HEALPix package. This falls rapidly as 
a function of I, as a result of beam smoothing. In order to correct 
for the smoothing, we divide the spectrum by the window func- 
tion of a symmetric Gaussian beam with FWHM=32'. This cor- 
rects for the mean smoothing and allows to recover the correct 
spectrum up to / max =300. Above that, the result begins to deviate 
from the input spectrum due to beam ellipticity, which cannot be 
corrected for by application of a window function only. 

Deconvolution recovers the TT input spectrum well up to 
Z=800. At higher multipoles there is little information left in the 
signal, and the convergence properties of the code become poor. 
The weaker EE and TE spectra are recovered up to /=400. Due to 
temperature leakage, the binned map only gives the correct EE 
spectrum for the 10 first multipoles in EE, and the TE spectrum 
is useless. 

We did one run with parameters / max =600, £ m ax=8 (not 
shown), to see if it would help to recover the EE spectrum to 
even higher multipoles. The results, however, did not differ sig- 
nificantly from the / max =600, k mdx -6 case. 
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Signal 




Error spectrum, signal 



1000 



Fig. 7. Noise-free 30GHz TT (top), EE (middle), and TE 
(bottom) spectra. Shown are the input (black) and the spec- 
tra constructed from the axim coefficients from deconvolution. 
Deconvolution results are shown for Ai max =4 (purple) and £ max =6 
(red), and for / max =400, 500, 600, 700, 800. The best-case re- 
sults (7 max =800, £ max =6) are plotted with a thicker linetype. For 
comparison we show also the spectrum obtained from harmonic 
expansion of a naive binned map (green). The blue line (TT and 
EE) is obtained by dividing the latter by the window function 
of a symmetric Gaussian beam with FWHM=32'. The TE spec- 
tra have been binned over 10 multipoles to reduce scatter and to 
make the plot more readable. 
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Fig. 8. 30GHz error spectrum, for T (top) and E (bottom), with- 
out noise. See main text for definition. The deconvolution results 
are shown for the same combinations of parameters as in Fig. [7] 
The input spectrum is shown in black. Also shown is the error 
spectrum for a harmonic expansion of the binned map, corrected 
with a window function of a symmetric Gaussian beam (blue), 
or (T only) with an ideal window function (light blue). The ideal 
window function was defined as the ratio of the input spectrum 
and the spectrum of the uncorrected binned map. 



5.2.4. Reconstruction error of harmonic coefficients 

A comparison between the reconstructed spectrum and the input 
spectrum is not a fully reliable test of the correctness of the de- 
convolution result. It is possible to construct a set of harmonic 
coefficients which gives the correct spectrum, although the indi- 
vidual coefficients are incorrect. 

In order to assess directly the accuracy of the axi m coeffi- 
cients, we construct a quantity we call the error spectrum. It is 
computed as follows. We take the difference between the axim 
coefficients we get from deconvolution, and the corresponding 
input coefficients . As before, X stands for T or E. We com- 
pute the spectrum of this difference as 

! l 

c *' r = 27TT ^ ^° xlm ~ a ^*^ axlm ~ a xim) 

m=—l 

The resulting spectrum is a measure of the error in the recon- 
struction of the axim coefficients, as a function of scale. 

We show the error spectra of the 30GHz runs for T and E 
coefficients in Fig. [8] For comparison, we show the input TT or 
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EE spectrum in the same figure. The axi, n coefficients are recov- 
ered with good accuracy when the error spectrum is well below 
the input spectrum. The meaning of the blue line is the same 
as in Fig. [7] that is, we compute the harmonic expansion of the 
binned map and divide them by the window function of a sym- 
metric Gaussian beam (FWHM=32')- 

We made another comparison with an ideal window func- 
tion, which we computed as the ratio of the uncorrected spec- 
trum of the binned map, and the input spectrum (green and black 
lines in the top panel of Fig. [7]). The error spectrum for this is 
shown in light blue in Fig.[8](TT spectrum only). By construc- 
tion, applying the ideal window function to the spectrum of the 
binned map returns exactly the correct input spectrum. The er- 
ror spectrum, however, is not improved with respect to the result 
obtained with Gaussian window function, which shows that the 
failure of the binning procedure to recover the correct spectrum 
above / = 300 in Fig. [7] was not just a result of a badly chosen 
window function. 

The error spectrum reveals that the deconvolved dxim coef- 
ficients are more accurate than those obtained by harmonic ex- 
pansion of the binned map already at the lowest multipoles. 



T d(!COi:"rl" 1 



5.3. Simulations with noise 

The first simulation was quite idealised, since the data was noise- 
free. We made another simulation, where we added white noise 
on top of the TOD. In other aspects the simulation was identical 
to the first one. No correlated noise was added. We are implic- 
itly assuming that correlated noise components, if present, were 
removed to sufficient accuracy prior to deconvolution. 

As expected, noise made the deconvolution process more 
difficult, and the code did not converge for the highest val- 
ues of Z max . We reached full convergence for the combinations 
'max=600, &max=6, and Z max =700, £ max =4. Case Z max =700, k max =6 
did not converge anymore. 

Fig. [9] shows a zoom into a point source, the same as in 
Fig. |U in presence of white noise. The deconvolution param- 
eters were / max =600, £ max =6. The deconvolved map smoothed 
with the fiducial FWHM=32' is contaminated by distorted noise. 
Also, we see some ringing around the point source due to the 
moderately low Z max . We therefore smoothed the map further to 
FWHM=45', which was enough to remove the distortion. For 
comparison we also smoothed the binned map to comparable 
resolution, by applying a smoothing by another FWHM=32' to 
it. This, combined with the original beam, gives a total resolution 
of roughly FWHM=45'. These two maps are shown in the bot- 
tom row. Even after the additional smoothing, the point source in 
the binned map is elongated. Deconvolution restores the circular 
shape of the source. 

Figure [10] shows a zoom into the Q polarisation map. We 
show the same patch of the sky as in Fig. [6] The raw binned map 
(not shown) is fully noise-dominated. We show maps smoothed 
to FWHM=45'. Again, deconvolution removes the galactic leak- 
age that is apparent in the binned map. 

Figure Q~T] shows the spectra constructed from the decon- 
volved axim coefficients, for the parameter combinations for 
which the code converged fully. Again we show also the spec- 
trum of the binned map, and the same corrected for a symmetric 
Gaussian beam with FWHM=32'. Deconvolution recovers the 
true TT spectrum up to nearly I = 600, after which noise begins 
to dominate. The spectrum of the binned map begins to deviate 
already below I = 400. In the case of the EE spectrum, both 
deconvolution and binning are unable to recover the true spectra 
except for the very lowest multipoles. In the case of the TE spec- 



-0.0010 ^^^^K 0.0010 -0.0010 0.0010 

Fig. 9. A zoom into the 30GHz temperature map in presence 
of noise: binned (left), and deconvolved (right). In the top row 
we show the unsmoothed binned map, and the deconvolved map 
with FWHM=32' smoothing. The deconvolved map is distorted 
by ringing and amplification of noise. In the bottom row we show 
the same maps smoothed further to FWHM=45'. The additional 
smoothing removes the distortion. 



Q binned Q deconvolved 




Fig. 10. A zoom into the 30GHz Q polarisation map: binned 
(left) and deconvolved (right) in presence of noise. The maps 
were smoothed to FWHM=45'. As in the noiseless case, the 
smoothed deconvolved map does not exhibit the spurious leak- 
age signal in the galactic plane. 
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trum, deconvolution performs considerably better than binning, 
and the result follows the input spectrum up to I — 300, although 
it is noisy. 

Finally, in Fig. Q~2]we show the error spectra in presence of 
noise. Deconvolution recovers the ajim coefficients well up to 
I = 600. In the case of «£/„, coefficients, the error is well above 
the input spectrum for all but the lowest multipoles. 

It is a well known and often quoted fact that deconvolution 
tends to amplify noise at high multipoles. We see this clearly in 
the upper deconvolved map in Fig. [9] The effects can, however, 
be suppressed by applying a sufficiently strong smoothing to the 
map, as we see from the lower panels of the same figure. 

We can look at the amplification of noise also in harmonic 
space f Fig.fTTb. While the TT spectrum of the binned map (green 
line) continues to fall until it levels out at I = 900 (not seen in the 
plot), the deconvolution spectrum begins to rise above / = 600. 
Note, however, that for the whole multipole range where we have 
a deconvolution result, the extra power in the deconvolved map 
is still below that in the binned map corrected for a symmetric 
beam (blue line). 

In the EE spectrum, the power of the deconvolved map ex- 
ceeds that in the binned map corrected for a symmetric beam, 
around I = 300. This is in the domain where we have full noise 
domination in any case. 

5.4. 70GHz simulation 

In order to test the performance of the code at higher multi- 
poles, we generated one noise-free simulation mimicking the 
data voulme of the Planck 70GHz channel. The used elliptic 
beams have an average FWHM of 13'; for more parameters 
see Table Q] We ran the deconvolution code with parameters 
'max =1700, k m . dx =6. 

Temperature and polarisation maps constructed from the 
axim with resolution FWHM=13' (T), or FWHM=18' (Q) are 
shown in Figs. [13] and [14] together with the corresponding 
binned maps. We show the same region of the sky as with 30GHz 
simulations. Again, deconvolution has worked well, returning 
the circular shape of the point source, and removing the galactic 
leakage in the Q polarisation map. 

In Fig. [15] we show the TT and EE error spectra for the 
70GHz simulation. The ajim are recovered well nearly up to 
I - 1700, and the a^im coefficients well above / = 1000. Thus 
we have shown that, with sufficient computational resources, our 
deconvolution method can be applied to realistic data sets of vol- 
ume and resolution of Planck 70GHz channel. 

5.5. Convergence 

The convergence criterion we had set in our code is rather strict. 
We require that the square of the residual vector has fallen below 
a fraction of 10~ 12 of its initial value. 

To see if a less strict criterion would help in reducing the 
computation time, we studied the convergence of the solution. 
We dumped the axi m coefficients after every 50 steps, and com- 
puted the error spectrum with respect to the input coefficients. 

We show the TT and EE spectra after 200 and 500 iteration 
steps in Fig. [15] along with the final results. We observe that the 
solution changes only very little after the first few hundred itera- 
tions, and 200 steps already give a significant improvement over 
the harmonic expansion of the binned map. In fact, after the first 
100 iteration steps we can see no visible improvement in the re- 
constructed sky map. We could thus reduce the computation time 
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Fig. 11. 30GHz TT, EE, and TE spectra, with white noise 
included. The TE spectra have been binned over 20 multi- 
poles. The deconvolution results are shown for three combi- 
nations of parameters: £ max =4 and Z max =600,700 (purple), and 
&max=6,/max=6 (red). The spectrum of a naive binned map is 
shown in green. The blue line shows the spectrum of the binned 
map corrected for a symmetric beam. The corresponding noise- 
free result is shown with dashed linetype in the TT plot, to indi- 
cate the region where noise begins to dominate over signal. 
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Fig. 12. 30GHz error spectra, in presence of white noise. The 
purple and red lines correspond to the same combinations of 
deconvolution parameters as inQT| The blue line shows the error 
spectrum for the harmonic expansion of a binned map, corrected 
for a symmetric beam with FWHM=32'. 



T deconvolved 





Fig. 13. 70GHz temperature map, without noise: binned (left) 
and deconvolved (right). The deconvolved map was smoothed 




Fig. 14. 70GHz Q polarisation map, without noise: binned (left) 
and deconvolved (right). Both maps were smoothed to final res- 
olution FWHM=18'. Shown is the same patch of sky as in Fig. [6] 
Deconvolution removes the spurious galactic signal. 



70 GHz error spectrum 




Fig. 15. 70GHz TT (top) and EE (bottom) error spectra, with- 
out noise. The blue line represents harmonic expansion of the 
binned map, corrected for a symmetric beam with FWHM = 13'. 



to resolution FWHM=13', corresponding to the mean width of The fully converged deconvolution result (/ m ax=1700, &max=6), 



the actual beam. Shown is the same patch of sky as in Fig. [4] 



is shown by a red line. We show also the error after 200 (purple 
dashed) and 500 (purple solid) iteration steps. 
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Fig. 16. Unconventional beam pattern 



binned deconvolved 




Fig. 17. Sky seen through a D-beam: binned (left) and decon- 
volved (fight) temperature map. The beam transforms the image 
of a point source into an image of the beam. Deconvolution re- 
covers the original shapes of the sources. The size of the shown 
patch is 83.3° squared. The horizontal structure near the bottom 
of the plots is the accumulation of weak point sources in the 
galactic plane. 



by a large factor, by using a more relaxed convergence criterion, 
without significant loss of accuracy. 

5.6. Complicated beam shape 

In order to demonstrate the power of the deconvolution method, 
we made yet another simulation with an unconventional D- 
shaped beam (D for "Deconvolution"), which is shown in 
Fig. [16] The linear dimension of the beam was approximately 
12° squared. For demonstration purposes we included only point 
sources and considered temperature only. For the binning pro- 
cess into the 3D map, we chose N^de = 1024 and = 256. 
Deconvolution parameters were / max = 500, k max = 30, and only 
a single full-sky survey was simulated. The comparatively high 
value of £ max was required in order to properly take into account 
the complicated azimuthal structure of the beam. 

A map binned directly from the TOD is shown in the left 
panel of Fig. [17] The beam transforms the image of a point 
source into a D-shaped pattern, resembling the beam. The decon- 
volved map is shown on the right; it was smoothed by a symmet- 
ric Gaussian beam with FWHM=1° in order to reduce ringing. 
Even with this complicated beam shape we are able to recover 
the underlying point sources. 

It is interesting that in this case it was sufficient to smooth 
the map with a beam of FWHM=1° to get a very satisfactory 
map, although the full dimension of the beam is much larger. 



6. Conclusions 

We have presented an efficient beam deconvolution code, de- 
signed for absolute CMB experiments, and tested it on sim- 
ulated data. The code is released under the terms of the 
GNU General Public License and can be obtained from 
http://sourceforge.net/projects/art-deco/ (as soon 
as the paper has been accepted). 

We have looked at maps constructed from the recovered aj\ m , 
aEim, and asim coefficients, and examined the coefficients di- 
rectly in harmonic space. These reflect different aspects of the 
solution. We compared the results to a harmonic expansion of a 
binned map. We have also shown that with sufficient computa- 
tional resources, we can extend the method to Z max =1700, which 
would be sufficient for the Planck 70GHz channel. 

In absence of noise we could recover the aj\ m coefficients 
to a high accuracy, and remove the effects of beam asymmetry. 
When white noise was added, we were able to reach a lower 
value of / max , but the advantage of deconvolution over binning 
was still clear. 

In the case of polarisation, deconvolution worked well in ab- 
sence of noise. We were able to almost completely remove the 
temperature leakage due to beam mismatch. When noise was 
added, results were less clear. Deconvolution removed the visi- 
ble galactic residual that arises from beam mismatch, but in har- 
monic space deconvolution did not seem to bring a clear benefit 
over binning. The recovered asim coefficients were dominated by 
noise at nearly all multipoles, but this was true for both decon- 
volved and binned maps. 

Our simulation used beams which were highly asymmetrical 
and also had a strong mismatch between them, which leads to a 
significant polarisation leakage. We did this in order to show the 
beam effects more clearly, but at the same time we also made the 
deconvolution problem very challenging. The accuracy we can 
expect when the code is applied to a real experiment, strongly 
depends on the beam shapes of the particular experiment. 

Some topics for future study can be mentioned. We have set 
a very strict convergence criterion for the CG iteration. We did 
a convergence study which indicates that a much more relaxed 
criterion could be sufficient for practical purposes. A future im- 
provement would be to define a more suitable convergence cri- 
terion. This could reduce the computation time further by a sig- 
nificant factor. 

In this paper we applied the method to a simulation data set 
which provides full sky coverage. Our code is built on a for- 
malism which makes no assumptions on the sky coverage. In 
practice, good convergence requires nearly full coverage. This 
is intuitively evident, since we are solving the harmonic coeffi- 
cients which necessarily represent the complete celestial sphere. 
If the input data only covers a part of the sky, the coefficients 
cannot be well constrained. A straightforward extension of the 
method would be to insert a prior which constrains the solution 
in the region which is not covered by the data, but this is beyond 
the scope of this paper and is a subject for future study. 

Yet another topic for further study is the determination of the 
noise bias present in the TT and EE spectra at high multipoles. 

Though the simulations used in this work mimick some as- 
pects of the Planck mission, the parameters used do not reflect 
the properties of the actual instrument. The results presented in 
this paper are thus not representative of the sensitivity of the 
Planck experiment, and the authors do not represent the Planck 
collaboration in this context. 
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Appendix A: Definitions 

In this appendix we present definitions an d conventions used 
in the paper. We follow the conventions of Varshal ovich et alJ 
(1988). 



A.1. Harmonic expansion of the CMB radiation field 

The CMB temperature and polarisation fields can be presented 
by spin-weighted harmonic functions with spin s = 0, +2. Spin 
5 = represents the temperature field, while components s = +2 
represent polarisation. 

The temperature field on the celestial sphere is expanded as 



CO / 

/=() m=-l 



(A.1) 



The Q and U Stokes fields, which describe linear polarisation, 
can be written as 



CO / 



(Q ± iU)(0, cf,) = J] ^ a mm ±2 Y lm (A.2) 



1=0 m=-l 



Here s Yi m are spin-weighted harmonic functions, and the 
corresponding (complex) coefficients. The harmonic functions 
and coefficients obey the symmetry relation 



1. Rotate by angle if/ around the z-axis (in positive direction). 

2. Rotate by angle 6 around the initial y-axis. 

3. Rotate by angle <f> around the initial z-axis. 

The rotation is equivalent to rotating first by angle (f> around the 
z-axis, then by angle 9 around the new y-axis, and finally by an- 
gle if/ around the new z-axis. The harmonic coefficients in the 
rotated coordinate system are given by 



a'slm = X a sl'n' D m>m(<f>, ^ <A) 



(A.8) 



where D m , m are Wigner functions. 

A Wigner function may be split into a product of three terms, 
each of which only depends on one of the angles: 



(A.9) 



Functions d l ,(ff) are the reduced Wigner functions. They have 
the property 

d 1 mk {Q)=8 mk . (A. 10) 

They are real-valued and obey the symmetry relation 

4tW> = (-i) m+ V- m ,-*(0)- (A. 11) 

The symmetry relation for the full Wigner functions is 



(A. 12) 



Me,® =(-i) ffl - s _. 5 y / %„(0, 

a s im =(-l)'"~ s «*. s ,/,- m 



(A.3) A.3. Wigner expansion 



The harmonic functions are related to the Wigner functions 
D' sm (6,cf,,ifj) through 



s Y lm (6,<p)=F l D l *(<[>,0,Q), 



where 



Fi = 



21+1 



1/2 



(A.4) 



(A.5) 



It is convenient to define further 



{ailm + a-2lm) 



(A.6) 



O-Blm - - 7rX a 2\m ~ O-2/m) 
2l 

The CMB radiation field is expected to be a statistically 
isotropic Gaussian spin 0,2 field. The statistical properties of the 
harmonic expansion coefficients are then fully characterised by 
four spectra Cj, Cf, Cf, Cf, such that 

( a Tlm a Tl'm>) =C/ <5//' 
( a *Elm a EI'm') =Cf6n>6 mm > 

( a *BIm a Bl'm') =Cf6u>5 mm > (A. 7) 

( a Tlm a El'm') -C l 5ff 



A.2. Wigner functions and rotations on a sphere 

The harmonic coefficients transform under rotations as follows. 

Assume we are given the coefficients a,/ m in one coordinate 
system. We then rotate the coordinate system, as described by 
Euler angles 0, 0, ifr: 



Wigner functions D' mk (<p, 6, iff) for integer values of /, m, k are de- 
fined in domain V: 

V : 0<<p<2n, < 6 < n, 0<ifj<2n (A.13) 

Domain V can be interpreted as the space of Euler angles. The 
volume element in domain V is dV = sin 8d6d<pdifr, and V's total 
volume is 87T 2 . 

Any finite function f(cf>, 8, if/) defined in domain V may be 
expanded in a series of Wigner functions as 

CO / / 

m W = XZZ c ' mk D' mk (^ 0, if/) (A. 14) 

/=() m=-lk=-l 

where the coefficients c l mk are given by 



'ink 



21+1 



■2tt /vr r*2n 

srnOdO I dif/f((f>,e,iff)D l ^,e,ifi) 
o Jo 

(A.15) 

Expansion dA.14t can be understood as the three-dimensional 
equivalent of the spherical harmonics expansion. 

In particular, a product of two Wigner functions 
Dj* t (w)Z/ /it ,(w) can be expanded as a series of Wigner 
functions as 



d:j^)d'm 



(A. 16) 



An integral over a product of three Wigner functions is given by 
Wigner 3 j symbols as 



^(^') J D, / rt '(«') J D ) ', 2 1 : fo («')^' 



(A.17) 



=8tt 2 (-1) 



m'+k' 



I I' h\l I I' h 



m -m' ni2 j\ k —kf &2 
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(IVarshalovich et al.lll988l) . 

The following properties of the 3 j symbols help in restricting 
the range of summation. A 3j symbol is non-vanishing 



h h h 



*0 



k\ k 2 h 
only if the two conditions 

k\ + k 2 + h = 

and 

\h - h\ <h<h + h 



(A.18) 



(A. 19) 



(A.20) 



are fulfilled. We can thus eliminate the summation over k 2 and 
m 2 . Product dA. 16b becomes 



(A.21) 



x(-iy 



i'+k' 



11' l 2 \ I I' h 



m —m' m! — m)\k —k' k' — k 



A.4. Beam 

AAA . Detector signal and definition of beam coefficients 

Consider a detector measuring CMB temperature and linear po- 
larisation. 

We assume now, very generally, that the signal recorded by a 
detector, apart from noise, is proportional to the « ,/,„ coefficients. 
We define the beam coefficients b s ik as follows. When the detec- 
tor is at some (predefined) fiducial position and orientation (for 
instance at the north pole), the signal detected in the absence of 
noise is 

?fid = ^ a simb* dm - (A.22) 

It is clear from the definition that the beam coefficients must 
contain all information not only of the beam shape, but also of 
the possible polarisation sensitivity and orientation. 

We may then rotate the beam to another position and ori- 
entation, as described in Section IA.2I The rotation brings the 
beam onto location 8, <p in spherical coordinates, in an orienta- 
tion defined by angle if/. The harmonic coefficients transform un- 
der rotations according to ( IA.8b - The signal seen by the rotated 
detector is then 



t = Y J a si m K lk D l * k (cf,,eA)- 



(A.23) 



slink 



In the following we derive expressions for the beam coeffi- 
cients in some special cases. 

A.4.2. Non-polarised beam 

Consider first a detector with an ideal beam, beam width equal 
to zero, and no sensitivity to polarisation. 

The signal detected by the beam is directly given by ( IA. lb . 
We see readily, by comparing ( IA.U and ( IA.22t . that if the chosen 
fiducial beam position is at the north pole (0 = 0, (f> = 0), the 
beam coefficients must be 



Ideal non-polarised beam: 



bsim = s^L(0, 0)<5 s o = Fid A )5„, 



(A.24) 



A general non-polarised beam with sensitivity g(8, (p) gives 
for the 5 = components 



General non-polarised beam: 

b 0b n= ff g(8,<f>)oY; m (e,cf>) sin 9d6dcf> 



(A.25) 



=Fi jj g(8, 4>)D' mQ (4>,6,0) sin 8d8dcp 
and b s im = for s — ±2. 

A.4.3. Polarisation measurements 

A polarisation-sensitive detector requires careful treatment, be- 
cause the Stokes parameters Q and U are not defined at the pole. 
The beam coefficients can still be defined in a meaningful way, 
as we show in the following. 

Consider again a beam at a fiducial position and orientation 
at the north pole. A general beam may have different polarisation 
sensitivity at different locations. We therefore define the beam 
shape g(9, <f>) separately for temperature and polarisation. 

The signal detected may be written in a general form as 

f = jj sm0d8d(f>{gT(8,(t>)T(e,(f>)+ (A.26) 
+ g P (9, 4>) [Q(0, 0) cos(2,A(0, cf>)) + U(0, </>) sin(2tA(0, 0))] } 

Here if/(0,<p) is the local direction of polarisation sensitiv- 
ity, measured as an angle from the local meridian. Functions 
grid, (p), gp(0, <p), and if/(0, <f>) together define a general beam. 

We proceed to calculate the beam coefficients. We write the 
trigonometric functions in exponential form, to obtain 

t = JJ sin8d8d(f>{gT(8,(f>)T(8,(f>)+ (A.27) 

+ \gp{8, 0) [«2 " iU)QS, cf>)e i2 ^ + (Q + iU)(8, <p)e- a ^\ } 

We can now insert the harmonic expansions ( IA. It and ( IA.2b into 
( IA.27t and read the beam coefficients 

boik = jj gT(8,(f>) Y* k (8, </>) sin flrfftty 

b ±2lk = \ Jj g P (8, 4>) ±2 Y; k (8, fte^S S in ft/ftty (A.28) 

We then use definition ( IA.4I ) and write the harmonic functions 
in terms of Wigner functions. In terms of the reduced Wigner 
functions we have 



= F ie - ik *d l k _ s (8) 



(A.29) 



The i/'-term can be transferred inside the Wigner function. We 
may thus write 

General polarised beam: 

boik =Fi Jj g T (8, 0)£>[, o (0, 8, 0) sin 6d6d<p 

b ± m =y jj gp(8, <f>)D l k ^ 2 (cf>, 8, ifj(8, <p)) sin 8d8dcf> (A.30) 

The above formulas for the beam coefficients are completely 
general, as long as the beam response is linear. 

We now make the simplifying assumption that the polari- 
sation sensitivity is "in the same direction" everywhere on the 
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beam. By "same direction" we mean that the direction of polar- 
isation sensitivity at an arbitrary location is obtained by rotating 
the polarisation direction vector at the north pole to the desired 
location along a meridian. The local polarisation orientation an- 
gle then becomes 

iA(M) = <Apoi-0 (A.31) 

where if/ po i is the deviation of polarisation direction from (f> — 
at the north pole. 

At this point it is convenient to write the Wigner functions in 
terms of the reduced Wigner functions dA.9l ). The beam coeffi- 
cients become 

Beam with unique polarisation direction: 

boik=Fi J sm0d6d' ko (6) J g T (0,<p)e- ik *d<p (A.32) 

^=Y^ j sin ed9d' k T2 (0) j g P (e,<f>)e- iik±2 ^d<[> 

We readily see that for a symmetric beam the only non- 
vanishing coefficients are those with k + s — 0. Especially, for a 
perfect beam with zero width and perfect polarisation sensitivity 
(£p — St) we obtain 

Ideal polarised beam: 

boik = F,d' kfi (0) = F t 6 k ,o (A.33) 

b ±2lk = Y<*U(0)e ±<2 *- = ye* *"**.*! 

We have assumed normalisation J g T (Q,)dQ = 1. 

Appendix B: Methods used to accelerate 
computation 

B.I. Load balancing 

In an MPI-parallel code, good load balancing among the indi- 
vidual tasks is crucial. As was mentioned in Section 13.21 some 
quantities are distributed across tasks depending on their colat- 
itude 6, or the index m. Since the computational cost of the al- 
gorithm is not independent of 9 and m, but rather an unknown, 
smoothly varying function of these two numbers, it is advisable 
not to assign sequential ranges of those indices to individual MPI 
tasks, because choosing index ranges which distribute the load 
evenly is very hard in that scenario. 

To avoid this complication, we adopted a "round robin" strat- 
egy for data distribution. Assuming n MPI tasks numbered 
to n - 1, we assign to a particular task i the 8 indices 6j, 6i +n , 
6i+2n, ■ ■ ■ , and deal with the index m in analogous fashion. This 
results in near-perfect load balancing. 

B.2. Convolution optimisation 

A significant part of total CPU time is spent in the Fast Fourier 
transforms necessary for the convolution with the quantity N 
(see section l3~2l i. Fortunately, several measures can be taken to 
decrease the resource requirements. First of al l, we are making 
use o f the highly optimised FFTW0 library (Frigo &" John son! 
1 1 9981) . which provides close to optimal performance for FFTs 
of specific size. 

4 |http://fftvTorg1 



Furthermore, it is important to notice that the execution time 
for an FFT of length n depends sensitively on the prime fac- 
torisation of n; in general, an n composed of only small prime 
factors is much preferable. 

Since we are performing a convolution operation, we have 
some freedom in the choice of array dimensions; they can be 
increased and zero-padded if desired. In our case, the minimal 
array dimensions are 4/ max + 1 and 4fc max + 1, respectively; if nec- 
essary, we increase both of them to the next larger number which 
is a product of the prime factors 2, 3 and 5 only. Depending on 
the prime factorisation of the original array dimensions, this can 
decrease the run-time for this part of the code by integral factors. 

Finally, the quantities that are to be convolved exhibit the 
symmetry 

which allows us to employ the specialised "real-valued" FFT, 
resulting in another substantial speedup and a size reduction of 
the precomputed Fourier transform of N (see Section lB~5l l. 

B.3. Efficient computation of the reduced Wigner functions 

During code execution, the reduced Wigner functions d' mk (9) are 
required twice in every conjugate gradient iteration step. Storing 
them in memory is prohibitively expensive, so it is important 
to regenerate them efficiently whenever n eeded. For this pur- 
pose w e use code originally presented by iPrezeau & Reineckel 
(12010 1 which was extended to make use of the SSE2 instruc- 
tion set present in all modern Intel-compatible CPUs. We also 
make use of the symmetry properties of the d matrix to avoid 
redundant computations. 

B.4. Shortcuts for symmetrical beams 

If a beam used for deconvolution is symmetric with respect to a 
180° rotation around its axis (which is true for elliptical beams, 
for example), all its b„i m coefficients for odd m vanish by defini- 
tion. This fact is used to skip unnecessary computations as well 
as save space by not storing any array elements which are known 
to be zero. It also reduces the minimum second dimension of the 
convolution arrays from 4k max + 1 to 2k m . dx + 1, saving even more 
time on the FFT operations and reducing memory consumption. 

B.5. Optional precomputation 

The convolution step mentioned in section 13.21 requires the 
Fourier transform of array N, which is constant throughout a 
run. The code permits to precompute and store this quantity, 
which reduces total CPU time (because only two FFT opera- 
tions are then needed for each convolution step, instead of three). 
However, especially for runs with high Z max and k max this in- 
creases the total memory consumption considerably. So for sit- 
uations where the available main memory is insufficient to store 
this precomputed array, we provide the option of recomputing 
the Fourier-transformed N whenever necessary, approximately 
halving the memory consumption of the code, but also slowing 
it down by 5-15 percent. For the 70GHz runs described in this 
paper, we had to make use of this space-saving feature. 

B.6. Supplying an initial guess 

While not strictly a performance optimisation, the code allows 
the user to supply an initial guess from which to start the CG 
iteration, which is then used instead of a vector of zeroes. This 
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Table C.l. Resource usage for various deconvolution runs 





/ 

'max 


k 

•MTiax 


1 wall 


T 

A setup 


T 

1 iter 


A/- 




30a 


400 


6 


10:07 


1:13 


3.9 


140 


6.6 


30b 


500 


6 


25:49 


1:20 


5.5 


270 


8.1 


30c 


600 


4 


1:24:57 


1:07 


5.1 


995 


6.4 


30d 


600 


6 


1:08:53 


1:35 


6.9 


591 


9.6 


30e 


600 


8 


1:40:40 


2:24 


10.3 


573 


11.9 


30f 


700 


6 


3:41:13 


1:48 


8.5 


1567 


11.3 


30g 


800 


6 


5:10:41 


1:56 


10.6 


1759 


12.8 


D 


500 


30 


10:31:38 


8:34 


20.3 


1842 


17.6 


70a 


1700 


6 


63:51:56 


15:27 


78.3 


2923 


28.8 


70x 


1700 


6 


9:43:05 


5:58 


11.8 


2955 




70y 


1700 


6 


5:25:32 


7:20 


6.6 


2922 





Notes. r wa u denotes the total wall clock time of the run, ^etyp is the 
time required for the code startup (i.e. data input, computing the RHS, 
Wigner transform of the hit count and preconditioner), ri ter gives the 
average wall clock time for a single CG step in seconds, N iteT is the 
required number of iterations until convergence, and "Mem" is the total 
memory usage in GB. Please note the additional comments in the text. 

is particularly helpful if one tries to compute series of deconvo- 
lutions on the same input data set, but with different parameters. 
For example, if a solution for a deconvolution problem with a 
given Z max and fc max has already been obtained, it can be used as a 
good starting point for other calculations with higher Z max and/or 
^max- In any case, the convergence criterion will still be evaluated 
with respect to an all-zero a s i m , not the supplied guess. 

Appendix C: Resource usage 

Most computations discussed in this paper were performed on a 
computer with 64GB of RAM and 24 Intel Xeon E7450 cores. 
Since this computer is used interactively, we only ran our calcu- 
lations on 16 of these cores and ensured that interactive usage did 
not exceed the capability of the remaining 8 cores. Nevertheless 
the measured timings contain small uncertainties and should be 
interpreted as upper limits. 

Table IC.fl gives an overview over the parameters of the var- 
ious runs and their impact on resource usage. Runs 30a-g are 
the 30GHz calculations discussed in Sect. 15.21 runs 70a,x,y are 
related to the 70GHz high-resolution study (Sect. 15.41 ). and run 
"D" is the deconvolution with the D-shaped beam presented in 
Sect. 15. 61 For all runs except the D-beam, the beam was assumed 
to be symmetrical with respect to a 180° rotation around its axis. 
Also, for all runs except the 70GHz ones the FFT of array N was 
cached (see Section lB~5l ). 

For the 70GHz run on the test machine described above (run 
70a), each CG iteration step needed significantly more time than 
for the calculations using the 30GHz data set; in combination 
with the high number of iterations, this leads to an inconve- 
niently long computation time. As a consequence, deconvolu- 
tion runs of this magnitude are better suited for execution on 
massively parallel computers. We therefore repeated the com- 
putation on the Louhi supercomputer of CSC (Cray XT4/XT5) 
on 128 and 512 processor cores, respectively. The corresponding 
performance figures are listed as runs 70x and 70y in Table IC.fl 
memory consumption was not measured for these runs. 

The computational power of a single computational core is 
roughly the same on both computers used for our tests. Looking 



at the total wall clock time for the 70GHz runs, it becomes ob- 
vious that the scaling from 16 to 128 cores is close to ideal, but 
becomes much worse when going from 128 to 512 cores. In this 
range, MPI communication and redundant computations begin 
to dominate. For even larger deconvolution problems, we expect 
that this transition will occur at a higher number of MPI tasks. 

The numbers of iterations required for the 70GHz runs are 
not identical, as one might intuitively expect. This is caused by 
the fact that in an MPI-parallel code some computation steps - 
most notably summations and other reductions of quantities re- 
siding on all tasks - are not carried out in a strictly defined order, 
but depend on the relative timing of the tasks during such an op- 
eration, which is nondeterministic. Since floating-point algebra 
with finite precision is non-associative, the results of such oper- 
ations can differ slightly from run to run, causing the conjugate 
gradient solver to take different paths towards the solution, but 
of course ultimately solving the same numerical problem. 

Compared to earlier publications on deconvolution map- 
making, it is evident that both memory and CPU requirements 
of the presented algorithm are quite modest and allow deconvo- 
lution up to values of / max that were prohibitively expensive up 
to now. Making use of massively parallel computing clusters, we 
expect that the algorithm can be applied to data sets with even 
higher resolution than those presented here. 
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